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Abstract. The fast marching algorithm computes an approximate solution to 
the eikonal equation in 0{N log N) time, where the factor log A' is due to the 
administration of a priority queue. Recently, [YBS06) have suggested to use 
an untidy priority queue, reducing the overall complexity to 0{N) at the price 
of a small error in the computed solution. In this paper, we give an explicit 
estimate of the error introduced, which is based on a discrete comparison 
principle. This estimates implies in particular that the choice of an accuracy 
level that is independent of the speed function F results in the complexity 
bound C'(i^max/^min-'V). A numerical experiment illustrates this robustness 
problem for large ratios -Fmax/i^min- 



1. Introduction 

The fast marching method, introduced in |Set96| . is an efficient algorithm to 
compute the discrete sohition to the eikonal equation || VT(a:)|j F{x) = 1, with speed 
function F. The algorithm is closely related to Dijkstra's single-source shortest path 
algorithm, and computes the grid-solution, starting with the grid-points adjacent 
to the boundary, by traversing the computational domain along increasing values 
of T. The total complexity of this method is 0{N log N), where N denotes the 
number of grid-points. Here, the factor logiV comes from the administration of a 
priority queue. 

In |Tsi95j . Tsitsiklis suggested independently a similar algorithm. Other than 
[Set96| . the discretization of the eikonal equation in |Tsi95| was obtained in the 
framework of optimal control theory. Moreover, Tsitsiklis showed that a bucket 
sort technique, together with a slightly different discretization, may be used to 
compute an approximate solution to the eikonal equation in 0{N) time. 

Another approach has recently been suggested by |YBS06| : using an untidy 
priority queue in the fast marching method also reduces the complexity to optimal 
0{N). The idea is to use a bucket sort technique together with a quantization that 
does not distinguish between values of T within a small range. This way, while 
sorting the values of T becomes less expensive, an error is introduced, however, 
which should be of the same order as the local truncation error of the discretization. 
The central result of this paper is Theorem IH which gives a rigorous estimate of 
the error introduced. Before we come to the theorem, the discretization and some 
basic notation is introduced in Section [H Then, in Section [3 we briefly recall the 
fast marching method and the bucket sort technique, as proposed in |YBS06| . A 
discrete comparison principle and the error estimation are subject of Section HI 
Finally, an example is given, which illustrates the theoretical result. 



1991 Mathematics Subject Classification. 65N22,65N06,65N15. 

Key words and phrases, fast marching method, eikonal equation, distance function, untidy 
priority queue. 



2 



CHRISTIAN RASCH AND THOMAS SATZGER 



2. Discretization 
Given a closed subset F C [0, 1]^, we consider the eikonal equation 

||Vr(a;)||^^(x) = l (xe [0,l]'\r), T|r = 0, (1) 
with continuous speed function F : [0, 1]^ ^ M, such that there are real numbers 

Fmin, F^^- with 

< F„in F{x) s; F^ax Vx. (2) 
The computational domain [0,1]^ is endowed with a Cartesian mesh [ih,jh) 
(i, J = 0, . . . , n) of grid-spacing h — 1/n. We partition the index set {0, . . . , n}^ in 
the two disjoint sets Qd and To, where represents the discrete version of the 
set r. Finally, we use the scheme proposed in [RT92| to discretize ||Vr||, arriving 
at the discrete equation 

maxp-.-r,-i^+-T,0)^ + max(i^,7r,-i?+^r,0)2 = i^, {^,J)enn 

iO (3) 
Ti] = 0, e Yd- 

Here Dfi^T denotes the backward/forward finite-difference approximation of the 
partial derivative with respect to x, that is, 



^ - h ' ^ h 

where = T{ih,jh), Fij — F(ih,jh). Furthermore, we set Ti^n+i = = oo 

and Tn+i.j = T-ij = oo (i, j = 0, . . . , n) for notational convenience. 

3. The Algorithm 

The solution of ^ can be efficiently computed with the fast marching method 
(compare e.g. [Set99]), which is given as follows. First, all points in F^i are tagged 
as known. All points in flu that have a neighbor in F^i are tagged as trial. Based 
on the known values in F^i, trial values are computed for all trial points, using 
(O. We denote by J\f the set of trial points, the so-called narrow band. Then the 
following iteration is performed: 

(Al) Let {«, j) G M denote the point with the smallest trial value. 
(A2) Tag (i, j) as known, and remove it from J\f. 

(A3) Tag all neighbors of {i,j), which are not known, as trial (unless this hasn't 

already been done) and add them to M . 
(A4) Recompute the values of T for all trial neighbors of («,j), using ([3]) with 

the known values of T only. 
(A5) If A^T^ 0, goto (Al). 

As the value of T^- in ([3]) depends only on smaller values of T in the neighboring 
points, the fast marching method computes a solution to (j3|. 
In the narrow band A/", the following estimate holds. 

Lemma 1. In every step of the fast marching method, there holds: 

max Ti^ — min T,,- < 



-I- l'] LLLLLL -t 27 ^ 

{^,3)^^^ F„ 

Proof. Let Af denote the narrow band at the beginning of step (Al), and Af* the 
narrow band at the end of step (A4), after one cycle of the fast marching method. 
Let {i,j) G Af denote the point with the smallest trial value, and let {k,l) denote 
some trial neighbor of {i,j), with its trial value Tki computed in step (A4). Then 

^ ^ fTki ~ Tij\^ . h 

-FT ^ Z ^ Tki !^ mm Tpq + —. 

Fli \ h J P,«)eAA Fkl 
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If we inductively assume the assertion to be true for JV, then, by the last inequality, 

h 



max T„„ ^ min T„„ r -= — 

(except for the neighbors of all grid- function values of points in J\f remain 

unchanged). On the other hand, the computation of a trial value, which de- 
pends on the new known value cannot yield a smaller result than Ty. Thus 
min(p^g)g_V Tpq min(p_,)gjv^. Tpq, and we are done. □ 

The crucial point in the fast marching method is step (Al). Typically, a heap 
based priority queue is used in order to keep track of the grid-point with the cur- 
rently smallest trial value Tij. Here we follow the suggestion of [YBS06| to use 
a so-called untidy priority queue within the fast marching method to reach linear 
run-time. The idea is to quantize the trial values contained in Af, and to pick any 
value, that quantizes to the minimal value, instead of the minimal value. 

The untidy priority queue is organized as a circular array of ub + 1 "buckets" , 
denoted by -Bq, ■ • ■ , , to store the narrow band Af. Those buckets are dynamical 
arrays, implemented as a singly-linked lists with FIFO ordering. Every bucket 
should contain trial points with trial values differing for at most 5, where 5 > is a 
small, appropriately chosen quantity. For the choice of the point with the smallest 
trial value in step (Al), we make no distinction between points contained in one 
bucket. In detail, the bucket Br contains all trial points {i,j) S Af with 

r = lT,j/6\ mod {riB + 1), where S — — . (4) 

^min • nB 

By the choice of S it is ensured that only grid-points with trial values within the 
range of (5 belong to a common bucket, that is \Tij — Tki\ < 5 for all {i,j),{k,l) S Br, 
as we have, similar to Lemma [T] 

max Tii — min < + 6 — S ■ (ub + i) 



max Tij 



min Tij 



< 1, 



when using the untidy priority queue. On the other hand, the modulus operation 
in ^ has the effect that buckets, which are emptied during the algorithm, are 
refilled in subsequent cycles of the fast marching method for memory efficiency 
reasons. During the computation we have to keep track of the bucket Bs, which 
holds the grid-points with the smallest trial values. The untidy queue supports two 
operations, namely 

• Insertion of some point To insert some point («, j), we compute r by 
Q and attach at the end of bucket Br- 

• Deletion of the point {i,j) with the approximately minimal trial value: The 
index s should hold the number of the bucket with the smallest trial value. If 
Bs is empty, we search cyclically for the next non-empty bucket B^i , passing 
cyclically from bucket to bucket, substituting s by [s + IJ mod (ns + 1). 
If all buckets are empty, the method terminates. Otherwise, we return the 
first element {i,j) in Bg, and remove it from B^. 

Since the recomputation of some trial value Tki in step (A4) can yield a smaller 
value, it might become necessary to move {k, I) to another bucket. In our approach 
we decided to insert a grid-point once more if it gets a new trial-value T^i such 
that the corresponding bucket number changes. Additionally, if we encounter some 
already known point in the delete operation, we simply skip the point and pick the 
next one. This approach may have the effect that a grid-point appears up to four 
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times in the untidy priority queueQ However, this way we avoid the costly removal 
of trial points from a bucket if they would have to be transferred to another one. 

The complexity of the insert operation is obviously 0{1). However, the compu- 
tational cost of getting the next minimal element in the queue may not be constant 
in a single step, as it incorporates the search for a non-empty bucket. Nevertheless, 
we obtain linear complexity in total, as the following lemma shows. 

Lemma 2. The complexity of the modified fast marching method with ub buckets 
is 0{N + n ■ ub), where N = \VId\- 

Proof. In every cycle (Al)-(A5) of the fast marching method, one grid-point be- 
comes known, and after N cycles the method terminates. Steps (A3) and (A4) 
require 0(1) operations. But possibly several buckets have to be passed in the 
delete operation in step (A2). However, the solution of ([3]) is bounded from above 
by AI = Let us imagine for simplicity, that we would use a linear array 

Bq, Bi, B2, ■ ■ ■ of buckets for our untidy priority queue, disregarding the modulus 
operation in Then we could cope with [M/(5J buckets. Hence also with the 
circular array, at most M/S — 2/h- ub — 0{n ■ ns) buckets have to be traversed 
in the delete operation during the whole algorithm. □ 

4. Estimates 

Let us write in short Nij{T) for the approximation of \\VT{ih, jh)\\ as given in 
formula ([3]). Then the following comparison principle holds. 

Lemma 3. Let Sij and Tij denote two grid-functions with Nij{S)Fij ^ 1 and 
Nij{T)Fij ^ 1 for all {i,j) G i^D, respectively. Then 

max (5'.y-Ty)+s$ max {Su-TijY 

Proof. Assume that for some inner grid-point {i,j) G il^ we have Sij — Tij = 
max(ij) (S'ij — Tij) = p > 0, and that additionally (i,j) is the grid-point with 
the minimal value of Tij among all grid-points in fix? realizing that maximum. The 
solution X of the discrete equation 

max(X-T,_i.„X-r,+ij,0)'+max(X-T,j_i,X-T,,j+i,0)2 = — . (5) 

at {i,i) can be written as0 

min js • min(Ti_i,,, Ti+i,,) + t • min(Ti j_i,Ti.,+i) + -^ys^^^l . 

s+t=l,s,t>0 [ V J J' ^ 'J J ' p.. J 

Hence, by the assumption Nij{T)Fij ^ 1, we have Ty ^ X, and there are numbers 
s, t ^ with s -\-t = 1, and fJ.,!^ (z { — 1, 1}, such that 



Tij > s • r,;+^j + t ■ Ti^j+^ + — VsM^. (6) 

J^ij 

Similarly we deduce by Nij{S)Fij < 1 that 

Sij ^ s ■ Si+^^j + 1 ■ Sij+„ + — + f^- 

i'ij 

Taking the difference of the last two equations yields 

P = ^ij ~ "^ij ^ * ' i^i+t^J ~ -^i+Mj) + ^ ' i^ij+f ^ Ti^jj^y) ^ p. 



^ In general, we observe that every point is inserted two times on average. 

^ In [BR06) . a variational discretization of the eikonal equation directly leads to this formula. 
As easily verified by solving a scalar extreme value problem, X fulfills l5l l. 
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By the choiceofp, we would have S'i+^j—Ti+^j = piis > 0, and Sij+,j—Tij+^ = p 
at > 0. Together with ^ we deduce that has a neighbor Tki with Tki < Tij 
and Ski — Tki = P- By our assumption, Tki has to be a boundary point, thus the 
lemma is proved. □ 

With the help of the discrete comparison principle, we are now able to prove 
an estimate on the error caused by the inexact minimization in the untidy priority 
queue. 

Theorem 4. Let T denote the solution, computed with the modified fast marching 
method, and let T denote the exact solution of (0). Then, for all {i,j) G flo, 

\T,j\ h Fmin UB 

where + 1 denotes the number of buckets, and 6 — — the increment. 

Proof. We analyze the error due to the untidy priority queue. When some value 
Tij becomes known., then there might be a smaller value Tki in the same bucket 
that would have been accepted before Ty if an exact priority queue had been used, 
with Tki ^ Tij < Tki + S. Let X denote the solution of 

max(X - r,_i.,, X - T,+i,,, 0)2 + max(X - r,.,_i, X - T,,,+i, 0)^ = 

Of course, we have X ^ Tij. Since, as we have seen, all values Tki that are accepted 
after fulfill Tki > Tij — S, we also have X > — S. Hence, 

/i2 .... .... 



2 ^ max(Ty - Ti_i J , T,j - T^+i j , 0) + max(rij - T.j^i, Tij - Ti^j+i, 0) 



< [nmx{X - f,^i,j,X - f,+i.,, 0) + 6]^+ [max(X - Aj-i,X - T,.j+i, 0) + 

Thus Nij{f)Fij ^ 1, and as Nij{T)F,j = 1, we obtain from Lemma[3]that ^ f,j 
for all {i,j). On the other hand, we have 9 ■ Nij{f)F,j = N,j{e ■ f)Fij 1 with 
9 = {1 + A/2(5Fmax//i)~\ < 6* < 1. Once again. Lemma [3] yields 

^ Ty - ^ (1 - 0) • 4- + max {9fki - Tki)+ V2 ■ ■ 5 ■ f„. 

□ 

Remarks: 

(1) Of course, the above result can be generalized to arbitrary space dimen- 
sions. In the estimate on the relative error, the factor \/2 has then to be 
substituted by the square root of the space dimension. 

(2) As Lemma [5] and Theorem show, an increase in the number of buckets 
causes a lower error at the price of a higher complexity. Moreover, we 
expect a total error of at least 0{h) due to the discretization of the eikonal 
equation, and a complexity of at least 0{N), as N is the number of grid- 
points. Hence, the suitable choice for ns is ub = 0{n) — 0{l/h). A choice 
of 71b = ^max/^min ' n buckcts would yield an error bound independent of 
the speed ratio Fmax/^min, but a total complexity of ©(FmaxZ-Fniin • A^)- 
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(3) The untidy priority queue can be used to modify the fast marching method 
on triangulations [KS98] , the method for the solution of the eikonal equation 
on parametric manifolds |SK04j . or the ordered upwind method for more 
general Hamilton- Jacobi equations given in }SV01| . to reduce the run-time 
to 0{N). 



relative error vs. F /F 




10" 10' 10^ 10^ 10'' 10^ 10^ 



Figure 1. Relative error vs. the ratio Fmax / i^min for different numbers of 
buckets (from top to bottom we used 2, 8, 32, 128, 512 and 2048 buckets). 



running time vs. gridsize 




numlDer of grid points ,.10' 

Figure 2. Run-time compared to tfie original fast marching method. 

5. Examples 

The first example concerns the error estimation asserted in Theorem U] We 
used the speed function F{x) = (1 + + sin(27r |a;|) with r — Fmax/Fmin 
on = [0, 1]^ and the boundary value r(l,0) = 0. Figure [T] illustrates the linear 
dependence of the relative error {Tij — Tij) j Tij on the ratio F,nax/Fmin for different 
numbers of buckets. 

The second example compares the run-time of the algorithm with an untidy pri- 
ority queue to the original fast marching Method. We computed on O = [0, 1]^ 
with the speed function l/F sampled from a [0, 1]— uniform distributed random 
variable. Figure [2] shows an asymptotically linear run-time of the modified method, 
and a non-linear behavior of the original method. 
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